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ABSTRACT 

To improve radiative transfer calculations for inhomogeneous clouds, a consistent means of modeling inho- 
mogeneity is needed. One current method of modeling cloud inhomogeneity is through the use of fractal param- 
eters. This method is based on the supposition that cloud inhomogeneity over a large range of scales is related. 

An analysis technique named wavelet analysis provides a means of studying the multiscale nature of cloud 
inhomogeneity. In this paper, the authors discuss the analysis and modeling of cloud inhomogeneity through the 
use of wavelet analysis. 

Wavelet analysis as well as other windowed analysis techniques are used to study liquid water path (LWP) 
measurements obtained during the marine stratocumulus phase of the First ISCCP (International Satellite Cloud 
Climatology Project) Regional Experiment. Statistics obtained using analysis windows, which are translated to 
span the LWP dataset, are used to study the local (small scale) properties of the cloud field as well as their time 
dependence. The LWP data are transformed onto an orthogonal wavelet basis that represents the data as a number 
of times series. Each of these time series lies within a frequency band and has a mean frequency that is half the 
frequency of the previous band. Wavelet analysis combined with translated analysis windows reveals that the 
local standard deviation of each frequency band is correlated with the local standard deviation of the other 
frequency bands. The ratio between the standard deviation of adjacent frequency bands is 0.9 and remains 
constant with respect to time. This ratio defined as the variance coupling parameter is applicable to all of the 
frequency bands studied and appears to be related to the slope of the data’s power spectrum. 

Similar analyses are performed on two cloud inhomogeneity models, which use fractal-based concepts to 
introduce inhomogeneity into a uniform cloud field. The bounded cascade model does this by iteratively redis- 
tributing LWP at each scale using the value of the local mean. This model is reformulated into a wavelet 
multiresolution framework, thereby presenting a number of variants of the bounded cascade model. One variant 
introduced in this paper is the “variance coupled model,” which redistributes LWP using the local standard 
deviation and the variance coupling parameter. While the bounded cascade model provides an elegant two- 
parameter model for generating cloud inhomogeneity, the multiresolution framework provides more flexibility 
at the expense of model complexity. Comparisons are made with the results from the LWP data analysis to 
demonstrate both the strengths and weaknesses of these models. 


1. Introduction 

The effect of clouds on radiative transfer has been 
shown to have a significant influence on model results. 
The magnitude of the radiative forcing effect of clouds 
has been studied during the Earth Radiation Budget 
Experiment (ERBE) and the International Satellite 
Cloud Climatology Project (ISCCP). Comparisons be- 
tween satellite data and general circulation model re- 
sults have shown discrepancies that are due in part to 
how clouds are represented in these models. Currently, 
a cloud layer is modeled as a homogeneous plane-par- 
allel feature that occupies a certain percentage of a 
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model grid box. Although this provides a reasonable 
first-order approximation for the purpose of radiative 
transfer, significant errors can be introduced by not in- 
corporating cloud inhomogeneity. As shown by Harsh- 
vardhan and Randall (1985), cloud layers with iden- 
tical mean liquid water path (LWP) can have dramat- 
ically different albedos. Methods of calculating 
radiative transfer in inhomogeneous media have been 
developed by Davies (1978), Stephens (1988), Ga- 
briel et al. ( 1993), Evans ( 1993), and others and have 
been used to study the magnitude of bias due to cloud 
inhomogeneity (e.g., Cahalan 1989; Lovejoy et al. 
1990; Kobayashi 1991 ) . Data collected during the ma- 
rine stratocumulus phase of FIRE [First ISCCP (Inter- 
national Satellite Cloud Climatology Project) Regional 
Experiment] indicate that cloud inhomogeneity re- 
duces cloud albedo by 10% or more (Cahalan et al. 
1994a). Though observational errors in cloud albedo 
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often exceed 10%, such a bias can have a noticeable 
effect on climate model predictions (Cahalan and Wis- 
combe 1992) . Having an even larger impact on remote 
sensing, the effect of cloud inhomogeneity on radiation 
thus merits careful study. 

The purpose of this work is to model cloud inho- 
mogeneity based on LWP data derived from the sur- 
face-based microwave radiometer during the marine 
stratocumulus phase of FIRE (Snider 1988) with even- 
tual application to radiative transfer calculations. Al- 
though analysis of this LWP data has been done by 
Cahalan et al. (1994a), we extend their work by ap- 
plying the technique of wavelet analysis. This approach 
not only provides information about the amount of in- 
homogeneity present at different temporal and spatial 
scales but also provides a multiresolution framework 
for modeling cloud inhomogeneity. This wavelet mul- 
tiresolution framework is a generalization of Cahalan’ s 
bounded cascade (BC) model and provides a means of 
changing the statistical properties of the model yet 
maintaining the multiscale nature of the LWP inho- 
mogeneity. 

Recent approaches to the study of cloud inhomoge- 
neity have drawn upon the popular topic of fractals and 
nonlinear dynamics. The connection of cloud fields to 
fractal concepts was made by Lovejoy (1982), who 
observed that cloud fields have a constant perimeter 
fractal dimension of 1.35 over areas ranging from 1 to 
1000 km 2 . Hentschel and Procaccia (1984) linked the 
fractal dimension to the turbulent diffusion involved in 
cloud development. Within this framework it is as- 
sumed that certain properties are invariant with respect 
to scale and that it is possible to derive features at one 
scale knowing the properties of the cloud at a larger 
scale. The connection between scales does not provide 
a deterministic representation of the smaller scale but 
rather a possible realization of the smaller scale that is 
statistically consistent with observed large-scale vari- 
ations. This connection has also been observed in spec- 
tra of wind speed and potential temperature analyzed 
by Gage and Nastrom (1986), who also observed a 
change in scaling at a few hundred kilometers. The re- 
lationship between scales also changes for areas 
smaller than 1 km 2 (Cahalan and Joseph 1989), where 
the scaling exponents depend on cloud type and bright- 
ness threshold. 

Simple macrophysical cloud models have been de- 
veloped that implement scaling properties. These mod- 
els are phenomenological rather than based on the 
equations of fluid flow and thermodynamics and are 
designed to give a simple representation of the inho- 
mogeneity within the cloud for purposes of calculating 
the cloud’s radiative properties. Schertzer and Lovejoy 
( 1987 ) considered a general class of multiplicative cas- 
cades. Cahalan (1989) studied the radiative properties 
of a one-parameter multifractal generated by such a 
cascade. Davis et al. ( 1990) focused on the thick cloud 
limit of a similar scaling model having a single fractal 


dimension. Using LWP data obtained during the ma- 
rine stratocumulus phase of FIRE, Cahalan developed 
a two-parameter fractal model for studying inhomo- 
geneity effects (Cahalan and Snider 1990) . This model 
redistributes LWP within a uniform cloud through an 
iterative cascade process that generates inhomogeneity 
at all spatial scales. Here we reformulate this model 
within a multiresolution framework. This provides a 
broader framework for modeling cloud inhomogeneity 
and leads to a number of variants, such as the variance- 
coupled (VC) model introduced in this paper. These 
multiscale statistical models generate features at dif- 
ferent spatial scales and then combine these features 
using an inverse wavelet transform. 

To accurately model the cloud inhomogeneity pres- 
ent at different scales, we focus on correlations between 
features of different scales as a function of time and 
frequency (or the analogous function of position and 
wavenumber). We first present a conventional “win- 
dowed’ ’ analysis in which the time-dependent relation- 
ship is studied by calculating statistics for data within 
a window that is translated across the length of the 
entire dataset. The time dependence of the frequency 
component is studied by performing a Fourier trans- 
form on the data located within each analysis window. 
Although the windowed analysis provides useful in- 
formation about the temporal and frequency character- 
istics of the LWP dataset, it does not optimize the win- 
dow size to the frequency being analyzed. An opti- 
mized window function balances the localization of 
temporal information with the localization of frequency 
information. 

Wavelet analysis was designed to optimize the trade- 
off between temporal and frequency localization. This 
is done by defining an oscillatory function of finite 
length. Convolution of this function with the dataset 
provides a time-dependent analysis of the data over a 
given band of frequencies. By stretching or compress- 
ing this oscillatory function, the mean frequency of the 
frequency band extracted by the wavelet analysis is de- 
creased or increased, respectively. Windowed Fourier 
analysis can be used to extract similar information; 
however, the size of the analysis window is indepen- 
dent of the frequency being analyzed. Therefore, the 
low-frequency component, which may consist of only 
one oscillation within the sampling window, is not sam- 
pled to the same degree as the high-frequency com- 
ponent, which consists of many oscillations within the 
sampling window. For wavelet analysis the effective 
window size changes with the mean frequency of the 
band being analyzed so that a similar number of oscil- 
lations are sampled. The Daubechies wavelets used in 
this study (Daubechies 1988) have an additional com- 
putational advantage over windowed Fourier analysis 
in that an orthogonal decomposition of the dataset is 
made using an efficient iterative algorithm based on 
multiresolution analysis. This is not true of a windowed 
Fourier analysis since the efficiency of the fast Fourier 



15 August 1995 


GOLLMER ET AL. 


3015 


transform is reduced by including a tapered window to 
the analysis. Reconstruction of the original data is also 
difficult for a windowed Fourier analysis since over- 
lapping windows are needed to account for the infor- 
mation lost while using the tapered window. 

This paper is presented in five sections. Section 2 
provides a description of the windowed and wavelet 
analyses used in this work. Section 3 describes the 
LWP dataset and presents the results of the analyses on 
this dataset. In section 4 the BC model is recast within 
the multiresolution framework. Within this context the 
VC model is introduced and compared to the BC 
model. Finally, section 5 discusses the primary obser- 
vations made during the LWP analyses and their rela- 
tionship to the modeling of cloud inhomogeneity. 

2. Windowed analysis 

In this section we describe two approaches to deter- 
mine the inhomogeneity within different frequency 
bands as a function of time and frequency. Since the 
data are not statistically stationary, the global statistics 
of the dataset are not representative of the statistics at 
any particular time within the dataset. Therefore, it is 
necessary to select subsets of the data over which the 
stationarity assumption is approximately satisfied. In 
the first approach, the time dependence of the dataset 
is studied by translating a window of fixed width over 
the entire dataset and computing statistics as a function 
of the window location. We will find that these win- 
dowed statistics provide useful information about the 
time dependence of the LWP dataset. By performing a 
Fourier transform at each window location, information 
about the frequency content is also obtained. The sec- 
ond approach uses a different windowing technique 
named “wavelet analysis,” which partitions the data 
into a series of frequency bands. Since the wavelet 
transform retains time information, the time depen- 
dence of inhomogeneity within each frequency band 
can be studied. To facilitate comparison between the 
two approaches, statistics from the result of the wavelet 
transform are calculated using translated windows 
matching those of the first approach. This provides 
common ground for discussing the nature of inhomo- 
geneity present within different frequency ranges. 

a. Windowed statistics 

The windowed statistics calculated from the LWP 
dataset include the first three moments of the data: the 
mean, standard deviation, and skewness. Although 
these statistics do not completely determine the prob- 
ability density function (PDF) of the data within the 
analysis window, they do constrain the shape of the 
PDF and provide useful statistics for modeling inho- 
mogeneity. Since these statistics are calculated for a 
data window of limited width, they represent a local 
statistic rather than a global statistic of the dataset. Use 


of the term “mean,” for example, will imply a “local 
mean” unless otherwise noted. 

The local statistics calculated for the LWP dataset use 
64-min windows that are translated by 32-min incre- 
ments. Since the sampling rate of the LWP dataset is 1 
min, a window of this width provides a sufficient statis- 
tical sample. This window size also restricts the segment 
of analyzed data to a length that nearly satisfies the sta- 
tionarity assumption. Changes in the statistics of the LWP 
data are primarily due to a diurnal cycle, although a syn- 
optic trend is also present. It is assumed that the solar 
forcing of the cloud layer does not change significantly 
over the length of 1 h, and therefore, a 64-min window 
is small enough to ensure stationarity. The validity of the 
stationarity assumption is studied by adjusting the win- 
dow width from 0.5 h to 8 h, where the later width clearly 
violates the stationarity assumption. 

The local statistics are also studied using a composite 
PDF that uses information from all of the analysis win- 
dows. A composite PDF is used since binning noise for 
a 64-point windowed analysis may obscure trends 
within the PDF. A composite PDF is generated by stan- 
dardizing the data within each analysis window and 
then combining all of the standardized data points into 
a single PDF. A standardized data point W' is calcu- 
lated using the equation 


07 

where /x, and 07 are, respectively, the mean and stan- 
dard deviation of the ith analysis window and W is the 
original value of the data point. This composite PDF 
indicates how the local PDF compares to a normal dis- 
tribution. Since the mean is removed and the standard 
deviation is set to 1, only the skewness and higher- 
order moments affect the shape of this PDF. The use- 
fulness of the composite PDF is limited to the assump- 
tion that the higher-order moments do not significantly 
change in sign over the extent of the entire dataset. For 
example, if the data distribution has local positive 
skewness, the predominant weighting of the composite 
PDF is to the left of zero (left of the mean) . However, 
if a similar number of data windows having both pos- 
itively and negatively skewed statistics are combined, 
the composite PDF may be symmetric about zero, in- 
dicating a nonskewed distribution. The composite 
PDFs used in this study are generated using 200 equally 
spaced bins ranging from —4 to +4 standardized de- 
viations about the mean. 

Frequency information is calculated for each win- 
dow by performing a Fourier transform. To avoid prob- 
lems due to boundary conditions at the edges of the 
analysis window, a Parzen’s window is used. Neglect- 
ing to use a windowing function affects the relationship 
between frequencies since there is a false high-fre- 
quency component introduced into the analysis (Press 
et al. 1992). The plot of the power spectrum provides 
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Fig. 1 . Plot of the (a) D8 wavelet function, (b) its Fourier transform, (c) the corresponding scale function, 

and its (d) Fourier transform. 


a convenient means of studying the relationship be- 
tween energy and frequency. Since the plot of the 
power spectrum is log-log, the energy as well as the 
frequency is normalized to a reference. Therefore, the 
slope of this plot yields a dimensionless quantity. From 
turbulence theory a constant negative slope indicates 
that the energy contained in the large-scale features 
cascades into the smaller-scale features. At the smallest 
scale, the kinetic energy is dissipated away as heat. In 
our analysis we assume that the inhomogeneity present 
at each timescale is related to the energy of the power 
spectrum at its corresponding frequency range. Since 
the LWP data used in this analysis represent a vertically 
integrated quantity, this inhomogeneity represents the 
horizontal nonuniformity of such processes as entrain- 
ment and vertical motion. For modeling purposes the 
slope of the power spectrum provides an important pa- 
rameter for applying scaling properties. 

b. Wavelet analysis 

Wavelet analysis is used to separate features into dif- 
ferent frequency bands while maintaining temporal in- 
formation. This makes it possible not only to study the 


relationship between the relative amounts of energy 
present in each frequency band but also to study any 
temporal dependence that might exist in this relation- 
ship. This is done by performing an inner product be- 
tween the LWP data and each member of a wavelet 
family. Each member of a wavelet family is 

related to the others by a combination of translations n 
and scalings (stretching of the function) m. Since nei- 
ther of these operations affects the relative shape of the 
wavelet, this analysis extracts features that are similar 
in shape. This is well illustrated in Weng and Lau 
( 1994), where examples are given for both continuous 
and orthogonal wavelets. 

The wavelet family used for analysis in this paper is 
based on the D8 wavelet (Fig. la). This orthogonal 
wavelet family was developed by Daubechies (1988) 
and is defined by eight “taps” (Table 1). These taps 
are used as weighting coefficients to extract the highest- 
frequency content from the data. The transform 

7 

d\(t„n) = <V'u(L), s 0 (t„)) = X giSoitn + 0 

i=0 

(n is even) (2) 


Table 1 . Tap values g,- and h, for the Daubechies D8 wavelet and scale functions, respectively. 






1 






0 

1 

2 

3 

4 

5 

6 

7 

gi 

hi 

0.0106 

0.2304 

0.0329 

0.7148 



0.0280 

-0.1870 



0.2304 

-0.0106 


Taken from Daubechies (1992). 
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multiplies the eight taps g ,• to the data series s 0 (t„) in 
order to generate the highest frequency detail d x (t„ n ) 
of the data series. The value of n is restricted to even 
values since the family of orthogonal wavelet trans- 
forms are translated two data points with respect to 
each other. As a result, the series generated by this 
transform is only half the length of the original data. 
This transform analyzes only the data at the smallest 
scale, which corresponds to the highest-frequency 
range. When the transform in Eq. (2) is applied itera- 
tively to multiple scales of the data, analysis by the 
discrete values of the taps are equivalent to analysis 
done by the nearly continuous function seen in Fig. la 
(Daubechies 1988). Although the transform behaves 
as a high-pass filter at the smallest timescale, at any 
larger timescale the transform performs a bandpass fil- 
tering of the data as shown by the Fourier transform of 
the D8 wavelet in Fig. lb. 

In general, the taps for Daubechies wavelets are de- 
fined such that members of the wavelet family are or- 
thogonal with each other. This is true only if the func- 
tions are translated by increments of two data points 
and the scale is changed by a factor of 2. These restric- 
tions allow the data to be expressed as a series of dyadic 
scales (Weng and Lau 1994). Mahrt ( 1991 ) discusses 
the statistics related to analysis with the Haar function, 
which is the simplest orthogonal wavelet. As with the 
Haar function, the D8 wavelet preserves energy and 
therefore allows a comparison of statistics between 
scales. In this study the D8 wavelet is chosen over the 
Haar wavelet so that leakage between frequency bands 
is minimized. Since the D8 wavelet is not a perfect filter 
function, there still exists some frequency leakage rep- 
resented by the lobes at 3 and 5 Hz in Fig. lb. The 
improved frequency characteristic is at the expense of 
some loss of temporal localization. 

As seen in Fig. lc, the D8 scale function <f> mM (t„) is 
defined by reversing the order of the D8 wavelet taps 
and negating each even tap (Table 1). This function 
performs a weighted sum of the data points by the for- 
mula 

7 

S\ So (?„)) = X hjS 0 (t n + i ) 

(•=0 

(n is even), (3) 

which generates the smoothed data 5 1 (r, i/2 ). This data 
has a resolution that is reduced by a factor of 2, and 
therefore, the transform acts as a low-pass filter as 
shown by the Fourier transform of the D8 scale func- 
tion (Fig. Id). As with the orthogonal wavelet trans- 
form, the translations are incremented by 2, thereby 
generating a series that is half the length of the original 
dataset. 

Due to the orthogonal nature of the wavelet trans- 
form, any information lost by the wavelet transform is 
maintained by its associated scaling transform. There- 
fore, the two series generated by these two comple- 


mentary transforms can be used to perfectly reconstruct 
the original dataset through the use of the inverse trans- 
form 

4 

So(t„ + k) = X g2i+kd\(t„n ~ i ) + h 2i+ kS\(t n i2 ~ i) 

<=() 

(k = 0, 1) (n is even). (4) 

From this notation it is seen that the same taps used in 
the forward transform are used to determine the weigh- 
ing of particular values of the smoothed, and 

detail, d l (t n/2 ), data to perform the inverse transform. 
It is also evident that the even taps are used for recon- 
structing the even locations of the data series and the 
odd taps for the odd locations. 

Although a D8 wavelet transform can be defined 
for scales other than the smallest scale, multiresolu- 
tion analysis (Mallat 1989) combines the wavelet 
and scale transforms to iteratively generate the wave- 
let transforms for larger scales. Since the D8 scale 
transform decreases the number of data points by a 
factor of 2, this low-pass or averaging transform ef- 
fectively increases the scale of the resultant coeffi- 
cients by a factor of 2. Therefore, a wavelet transform 
performed on the resultant components of the scale 
transform s i (t rl , 2 ) is the same as performing the 
wavelet transform of the next largest scale on the 
original data s 0 ( t n ) . 

Notation for the multiresolution analysis with the D8 
wavelet is a generalization of Eqs. (2), (3), and (4). 
The scale transform 

7 

■S/nT 1 (^ 1 / 2 ) ( ( tn ) ■> ^ hjS m (t n "i“ 1 ) 

i=0 

( n is even) (5) 

reduces the resolution of the smoothed data s m ( t„ ) 
obtained from the previous level of the scale trans- 
form. Once again, notice that the length of the (m 
+ 1 )th-level scale transform is half the length of the 
mth-level scale transform. The rath-level scale trans- 
form spans the original data series with N/2 m points 
since translations are by increments of 2 and each 
change in scale increases the translation distance by 
an additional factor of 2. In this case N is the number 
of points present in the original data series. The 
wavelet transform 

7 

dm + li k t n / 2 ') ( tyltl ,H ( t,; ) , S m (t;;)) gi^lll ( hi “I" i) 

(=0 

(« is even) (6) 

removes the detail from the smoothed data obtained 
from the previous level of the scale transform. This 
iterative definition for the wavelet and scale trans- 
form makes computation very efficient since larger 
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Fig. 2. Diagram of multiresolution analysis. 


levels of the transform are performed on a reduced 
length of the dataset. The inverse function 

4 

+ k) = X 82 i+kd m + l(tn/2 ~ 0 
i—0 

+ h 2i + k s m+x {t,„ 2 -i) (k = 0,Un is even) (7) 

is also defined iteratively and provides perfect recon- 
struction given that infinite precision is maintained in 
the transform. 

This iterative transform process is illustrated in Fig. 
2, where the multiresolution analysis is performed on 
data sampled at 1-min intervals to generate smoothed 
and detail information at the 2-, 4-, and 8-min scales. 
The original LWP dataset, which has a resolution of 1 
min, is represented by the rectangle labeled .v 0 ( t „ ) . Per- 
forming the scale transform with the scaling function 
(f) l n reduces the resolution (increases the scale) of the 
data by a factor of 2 and generates the smoothed data 
s x (t n i 2 ). The index t nl2 is used since the transformed 



Fig. 3. Frequency ranges of multiresolution analysis 
for the LWP data. 


data have only half the number of points of the LWP 
data while spanning the same time period. Therefore, 
each data point in S\(t n , 2 ) represents a local 2-min 
average of the LWP dataset. Although an arbitrary 
position index of n could be used, Fig. 2 uses the index 
t Hl2 - to simplify the collocation of features from the 
scale m , with a location within the original LWP data- 
set. Since the wavelet and scale function complement 
each other, the detail lost by this reduction in resolution 
is retained by the detail d x (t n , 2 ) obtained through the 
wavelet transform. To obtain the next lowest resolution 
(next largest scale) the multiresolution process is re- 
peated on the transformed dataset S\(t„, 2 ) to obtain the 
smoothed data s 2 (t „, 4 ) and the detail data d 2 ( t „ /4 ) . At 
this scale features are representative of the 4-min time- 
scale and each transformed data series has only a quar- 
ter of the number of points of the LWP dataset. This 
process is repeated once again on s 2 (t n/4 ) to obtain both 
the smoothed, s 3 (r w/8 ), and the detail, d 2 (t n/& ), datasets 
at the 8-min timescale. 

Continuing this processes to the largest scale allows 
the original dataset to be represented by a series of 
wavelet transforms that span all possible scales. Since 
the scale transform can be represented by the scale and 
wavelet transform at the next largest scale, the concept 
of the scale transform serves only as an intermediate 
step for iteratively generating the wavelet transform at 
all scales. At the largest scale a single component rep- 
resents the scale transform, and its value corresponds 
to the global average of the LWP dataset. Since this 
analysis is based on an orthogonal decomposition, all 
information from the LWP dataset is retained, and it is 
possible to reconstruct the original dataset from the 
components of the wavelet transform at all scales. 

The multiresolution analysis represented in fre- 
quency space reveals the trade-off between frequency 
and temporal localization. Figure 3 illustrates the fre- 
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quency ranges represented by each scale of the wavelet 
analysis. At the smallest scale the wavelet transform 
acts as a high-pass filter where the cutoff frequency is 
at half of the Nyquist frequency. Since the sample rate 
of the LWP data is 1 min, the Nyquist frequency is 
30 h _l and the cutoff frequency of the first level detail 
d\(t n , 2 ) is 15 h '. The length of time represented by 
each component of the first level detail is 2 min. At the 
next largest scale each component of the wavelet anal- 
ysis represents 4 min, while in frequency space these 
components represent information that ranges between 
7.5 and 15 h -1 . In a time-frequency phase space rep- 
resentation the product of the time localization and the 
frequency localization remains constant. Therefore, as 
the scale of the wavelet analysis increases, lower fre- 
quencies are analyzed, the range of frequencies con- 
tributing to the analysis decreases, and the number of 
original data points contributing to the analysis in- 
creases. 

From this presentation of the multiresolution analy- 
sis an obvious limitation exists that affects data anal- 
ysis. Since the change in scale is limited to factors of 
2, a translation of the dataset by one position will 
change the representation of the wavelet transform of 
the data. In addition, if any one frequency has a strong 
influence on the dataset, the wavelet analysis will at- 
tribute it to a range of frequencies or, worse yet, to two 
frequency ranges when the frequency lies near a tran- 
sition between two scales. Since the scales analyzed by 
the wavelet transform are powers of 2 of the resolution 
of the original data, the analyzed scales cannot be 
shifted apart from remapping the LWP data into a dif- 
ferent sampling rate. Although many of these problems 
can be resolved by using nonorthogonal wavelets, these 
problems are not prohibitive to the analysis used in this 
research. Since statistics are calculated from the wave- 
let transform, any change in the specific representation 
does not significantly change the statistical results. The 
features present in the LWP data used in this research 
are the result of turbulent processes that do not favor a 
specific frequency over another. The only predominant 
frequency present in the data is the diurnal cycle, and 
any analysis at that scale violates the stationarity as- 
sumption. The advantage of retaining the orthogonal 
condition is that features present in the data cannot be 
overrepresented in the wavelet analysis. This provides 
grounds for comparing statistics contributed by differ- 
ent scales without a bias due to one frequency having 
a stronger contribution than originally present in the 
data. 

The correspondence between windowed analysis and 
wavelet analysis is illustrated in Fig. 4. The window T„ 
is chosen to be about 1-h wide (64 points) and is trans- 
lated by about 0.5 h (32 points) increments. Although 
there is overlap between adjacent windows, this permits 
a larger number of windowed statistics to be generated. 
Of comparable size to the data windows, the fourth- 
level D8 wavelets t// 4 ,„ are translated by only 16 points. 



Fig. 4. Correspondence between windowed analysis 
and wavelet analysis. 


In addition to the difference in translation, the weight- 
ing of the data is equal over the window T„, while 
transforms with the */r 4> „ wavelets have their largest 
weighting based on the middle half of the window. An- 
alyzing inhomogeneity of the data over a window T„ 
includes all scales smaller than the window size. How- 
ever, by using the , and t/r 3 „ wavelet functions, 

inhomogeneity is partitioned into a series of scales that 
lie in the 32-, 16-, and 8-min scales, respectively. Al- 
though the size of the windowed analysis can be ad- 
justed to these different sizes, it does not provide an 
orthogonal representation between scales and therefore 
does not provide a simple method for reconstructing 
the original dataset. 

c. Wavelet statistics 

The important statistics calculated from the wavelet 
analysis are the standard deviation and the Markov pa- 
rameter for each level of detail. The wavelet transform 
partitions the inhomogeneity of the data into a number 
of time series each associated with a different fre- 
quency band. The windowed analysis discussed in sec- 
tion 2a is applied to each of the wavelet-generated time 
series to calculate a standard deviation. This statistic is 
equivalent to the square root of Gamage and Hagel- 
berg’s (1993) wavelet variance. This scale-dependent 
standard deviation is calculated for the first-level detail 
d\ (t, l/2 ) through the third-level detail <i 3 (f„ /8 ). Since 
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each level of detail has half the number of points as the 
preceding level, the number of points used for calcu- 
lating the detail standard deviation is halved. For a 64- 
point window, the standard deviation of the first-level 
detail d i (t„, 2 ) is calculated using 32 points, the second 
level using 16 points, and the third level using 8 points. 
Due to the limited statistics for calculating the standard 
deviation of the third-level detail, analysis focuses on 
the first two levels and the relationship between them. 

The ratio of mth-level standard deviation to (m 
+ l)th-level standard deviation is related to the slope 
of the power spectrum. Since the wavelet transform 
separates data into different frequency ranges, it is pos- 
sible to relate a particular level of the detail to a par- 
ticular portion of the Fourier spectrum. For data that 
has a linear power spectrum, the ratio of standard de- 
viation between adjacent levels of the wavelet trans- 
form should be constant. This in turn is related to the 
self-similar or self-affine nature of the data. The rela- 
tionship between the ratio of standard deviation 
between levels and slope of the power spectrum is 
based on observations while conducting our research 
and still needs to be quantified by a mathematical re- 
lationship. 

An additional statistic for the mth-level detail is the 
Markov parameter. This parameter is equal to the lag 
1 autocorrelation of the mth-level detail. Since each 
level of the detail spans a limited frequency band, the 
Markov parameter measures the correlation of inho- 
mogeneity at a particular scale. The wavelet transform 
is not a perfect filter, and therefore, it is possible for 
structures at the next largest scale to generate autocor- 
relation at the level being analyzed. The presence of 
autocorrelation at all scales affects the slope of the 
power spectrum and complicates the mathematical re- 
lationship between the ratio of standard deviation be- 
tween levels and the slope of the power spectrum. The 
presence of negative autocorrelation at all scales causes 
a time series to be more intermittent than one with pos- 
itive autocorrelation given that the slope of the power 
spectrum is forced to be the same. The LWP dataset 
used in this analysis has intermittent features as re- 
ported by Davis et al. ( 1994) and, therefore, makes the 
Markov parameter important for modeling applica- 
tions. The term autocorrelation or Markov parameter is 
used interchangeably, since only the lag 1 autocorre- 
lation is used in this study. 

3. Vertically integrated liquid water data 

a. Description 

The data used for this study were obtained during 
the marine stratocumulus phase of FIRE. The intensive 
field operation was conducted from 1 July through 19 
July 1987 and consisted of ground-based, airborne, and 
satellite measurements. The ground-based instruments 
were stationed on San Nicolas Island, which is located 


100 km southwest of Los Angeles. One goal of the 
experiment was to improve the characterization of 
clouds for radiative transfer studies (Cox et al. 1987). 
The LWP data are well suited for this since they pro- 
vide a lengthy record of cloud inhomogeneity and con- 
tain changes therein due to time development. 

The LWP data used in this study were obtained from 
a ground-based, upward-looking microwave radiome- 
ter. Derived from an algorithm using the 20.6- and 
31.65-GHz channels, the LWP data have an absolute 
accuracy of 20% (Albrecht et al. 1990). Since the re- 
trieval method may contain a bias that affects the ab- 
solute amount of LWP, retrieval of variations in the 
LWP can be known to a greater accuracy of 10% (Ca- 
halan and Snider 1989). Consisting of 1-min averages 
of 10-Hz data, this LWP dataset runs nearly continu- 
ously for 19 days. Missing data gaps of several hours 
are present in the data along with 2-min gaps approx- 
imately every 4 h, which are used for calibration of the 
equipment. 

Since the LWP measurements are from a zenith view 
of the sky, the data series records not only time devel- 
opment during the course of the 19 days but also spatial 
development as the clouds advect over the instrument. 
The instantaneous sampled volume for the microwave 
radiometer is a truncated cone that has an average di- 
ameter of 21 m and depth of 100 m, assuming an av- 
erage cloud base of 425 m during the experimental pe- 
riod (Blaskovic et al. 1991). Assuming an average 
wind speed at cloud height of 5 m s" 1 , the 1-min mea- 
surements of the microwave radiometer are an average 
of LWP over a 300 m X 21 m X 100 m section of 
cloud. Using 64-min analysis windows, the spatial ex- 
tent for each window corresponds to a 19.2-km section 
of cloud. On this timescale the data are assumed to be 
stationary as discussed in section 2a, and therefore, 
variability within the window is assumed to be the re- 
sult of spatial structure rather than time development 
of the cloud field. 

b. Analysis 

The LWP data contain unusable portions due to 
missing data and clear skies. Using only data from 9 
through 17 July eliminates most of these unusable sec- 
tions. The bias due to the remaining missing data is 
reduced by linearly interpolating over these short gaps. 
Since the greatest variability is found on timescales 
long compared to the gaps, the interpolation introduces 
little error. Observations of LWP below a threshold of 
10 g m~ 2 are assumed to be unreliable due to the bias 
present in the microwave retrieval algorithm for thin 
clouds and due to the presence of clear skies. To avoid 
problems due to poor sampling, any window containing 
more than 25% clear skies, thin clouds, or missing data 
is rejected. When calculating the mean, standard de- 
viation, skewness, and PDF, only cloudy points are 
used within each window. However, when calculating 
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Fig. 5. Time development of the LWP data’s local (a) mean and (b) slope of the power spectrum. 


the Fourier transform and any statistic related to the 
wavelet transform, all data points are used since con- 
tinuity between adjacent points is important. For these 
last two calculations, the original values are used since 
it is expected that thresholding the data will introduce 
a false high-frequency component to the data. 

The mean of the LWP has a strong diurnal cycle as 
seen in Fig. 5a. The mean is greatest just before sunrise 
and reaches a minimum in the late afternoon. This vari- 
ation in LWP is related to changes in both the cloud 
base and cloud thickness (Albrecht et al. 1990) . A syn- 
optic trend is also present in the data as demonstrated 
by the differences in the daily maximums of the LWP. 

Although the mean of the LWP has a strong diurnal 
dependence, the slope of the power spectrum (Fig. 5b) 
is nearly time independent. As discussed in section 2a, 
the slope of the power spectrum is related to the fre- 
quency dependence of energy and therefore is related 


to the scale dependence of inhomogeneity. There are 
fluctuations in the value of the power spectral slope; 
however, there is no obvious trend in this statistic. Us- 
ing 1-h analysis windows the slope has a mean value 
of —1.75. This value remains at —1.75 for windows 
that extend to 8 h in length, although a window of this 
size violates the assumption of stationarity. For window 
sizes beyond 8 h, the slope of the power spectrum is 
affected by the diurnal cycle; however, the value of the 
slope fluctuates about the same global value of —1.75. 

The most important result of the LWP’s wavelet 
analysis is the relationship between inhomogeneity at 
different scales. The degree of inhomogeneity is related 
to the fluctuations in the LWP and is quantified through 
the statistic of standard deviation. The standard devia- 
tion for the LWP data (Fig. 6a) consists of a diurnal 
cycle with a daily peak value near 1200 UTC, which 
corresponds to dawn. The standard deviation decreases 
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Fig. 6. Time development of the LWP data’s local (a) standard deviation and standard deviation 
of the level (b) 1 , (c) 2, and (d) 3 wavelet transform. 
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during the day to reach a minimum in late afternoon. 
Due to the sampling rate of the microwave radiometer, 
the frequencies present in this dataset range from 0 to 
30 h~‘. The first three levels of the wavelet analysis 
partition the original dataset into the frequency ranges 
of 15-30, 7.5-15, and 3.75-7.5 h _1 . The standard de- 
viation of these three levels (Figs. 6b-d, respectively) 
corresponds well with the standard deviation of the 
original dataset and with each other. The windowed 
standard deviation also contains a synoptic feature that 
has a time period of approximately 7.5 days. This syn- 
optic feature also exists in the standard deviation of the 
wavelet transform for levels 1, 2, and 3. 

The correlation between inhomogeneity at different 
scales is very strong and is not influenced by the diurnal 
cycle. The correlation coefficient between level 1 and 
level 2 standard deviation is 0.81. This implies that 
66% of the changes in local standard deviation at the 
2-min scale is explained by the changes in local stan- 
dard deviation of the 4-min scale. At higher levels of 
the wavelet transform the correlation coefficient be- 
tween adjacent levels of standard deviation is also near 
0.81. The ratio between level 1 and level 2 standard 
deviation (Fig. 7a) is also time independent, although 
there may be a dependence on the synoptic scale as 
indicated by the slight increase in the ratio around days 
13 and 14. This lack of time dependence is also ob- 
served in the ratio of standard deviation between levels 
2 and 3. Another important feature is that the ratio of 
standard deviation is independent of which two adja- 
cent scales are being compared. 

The ratio of standard deviation between scales is ap- 
proximately 0.9 and is valid up to the ratio of standard 
deviation between levels 5 and 6. The ratio of standard 
deviation between levels 6 and 7 drops to 0.7 and is 
most likely due to the size of the analyzing window. 
Although the ratio of standard deviation between levels 
1 and 2 uses a 1-h window, the window size is doubled 
at each increase of scale to ensure sufficient statistics. 
The ratio of standard deviation between levels 6 and 7 
depends on a window size of 32 h, which is influenced 
by the synoptic trend of the data. Each wavelet coef- 
ficient within this window is representative of 1 h for 


level 6 and 2 h for level 7. Therefore, the reduction in 
the ratio of standard deviation between scales may be 
a combined effect of the diurnal and synoptic trends of 
the data. It is interesting that the ratio of standard de- 
viation between scales remains unaffected by window 
sizes up to 32 h, which is well beyond a reasonable 
stationarity assumption. 

As discussed previously, the windowed mean (Fig. 
5a) has a strong diurnal dependence, and this depen- 
dence is positively correlated with the windowed stan- 
dard deviation (Fig. 6a) by a correlation coefficient of 
0.67. However, on the 2-min scale the correlation co- 
efficient between the mean and the level 1 standard 
deviation is 0.56. This means that only 30% of the 
change in standard deviation at the 2-min scale is ex- 
plained by the mean. This is made evident by taking 
the ratio between the level 1 standard deviation and the 
mean (Fig. 7b). This ratio is affected by the diurnal 
cycle and has a minimum of 0.03 at dawn and a max- 
imum of 0. 1 in late afternoon. This ratio does not ap- 
pear to contain a synoptic trend and therefore is useful 
for modeling as long as the diurnal dependence is taken 
into account. Although the ratio of standard deviation 
between scales is a better predictor of small-scale in- 
homogeneity, the ratio of standard deviation to mean 
is useful when information about the inhomogeneity of 
the LWP is unavailable. 

The windowed Markov parameter for the first level 
of the wavelet transform does not contain a diurnal cy- 
cle. Having a mean value of —0.28, this statistic be- 
comes less negative around 13 July, when the mean 
LWP is small. This statistic is also time independent at 
the second and third levels of the wavelet transform 
and has mean values of —0.28 and —0.29, respectively. 
Beyond these levels the statistic becomes less negative 
and may be affected by diurnal influences. The auto- 
correlation at each level of the wavelet transform in- 
dicates that correlated structures exist at all the ana- 
lyzed scales. This observation corresponds well with 
the fractal concepts of self-similarity and self-affinity. 
Since wavelet analysis restricts the inhomogeneity to a 
given timescale, the autocorrelation is only indirectly 
related to larger-scale structures within the cloud field. 



July 87 



July 87 


Fig. 7. Time development of the LWP data’s local (a) ratio of standard deviation between levels 1 and 2, 
and (b) ratio between level 1 standard deviation and mean. 
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Due to the existence of the Markov parameter, it is 
important that any inhomogeneity model manifest this 
time-independent statistic. 

Unlike the statistics discussed previously, the skew- 
ness of the LWP is both time and scale dependent with- 
out any obvious connection to the time dependence of 
the mean. From the 1-h analysis windows the skewness 
fluctuates between — 1 and + 1 and has a mean skew- 
ness of 0. 1 . When the analysis window size increases 
to 2 h, the skewness increases to 0.28. Skewness at this 
scale is less uniform and appears to have some time 
dependence to it. Extending the analysis window to 4 
h gives a global skewness of 0.41 and begins to reveal 
the diurnal dependence of the skewness. The skewness 
over this larger window has a weak diurnal dependence 
that is also related to the mean LWP. When the mean 
LWP remains relatively constant, the skewness is zero 
or even negative. However, when the local LWP begins 
to change, the skewness is positive. This appears to be 
related to the nonlinear way in which cloud thickness 
develops. When the cloud LWP increases, it does so at 
a faster than linear rate. As a result, the PDF of LWP 
over a set interval develops a long, positively skewed 
tail. Likewise, as the cloud decreases in LWP, the dis- 
sipation is faster than linear and the distribution of 
LWP for a set interval is again weighted toward a thin- 
ner section of cloud with a few thick sections of cloud. 

The dependence of local skewness on analysis win- 
dow size is best illustrated using the composite PDFs 
for different window sizes. The composite PDF for a 
1-h window (Fig. 8a) appears nearly Gaussian when 
compared to a normal distribution. However, as the 
window size increases, the largest weighting of the 
PDF shifts to the left and begins to generate a longer 
positive tail. This shift in the PDF is the result of vio- 
lating the stationarity assumption since the statistics of 
the cloud change significantly over 8 h due to diurnal 
effects. For window sizes of 8 h (Fig. 8b) and larger, 
the PDF becomes more lognormal, as reported by Ca- 
halan and Snider ( 1990) . Although the contribution of 
positive skewness is present over the entire dataset, a 
large amount of skewness exists on 17 July during a 
drizzle event where the mean LWP becomes very large 
over a short period of time. 



Standardized Deviation 

Fig. 8. LWP composite PDFs 


Since it can be argued that the PDF of the 1-h win- 
dow is dominated by instrument noise and thereby pre- 
dominantly Gaussian, an alternate composite PDF is 
generated. The alternate PDF uses only windows that 
have a standard deviation at least four times as large as 
the expected standard deviation due only to noise. The 
skewness dependence of the 1-, 2-, 4-, and 8-h alternate 
PDFs is the same as the composite PDFs and therefore 
indicates that this observation is not due to noise con- 
tamination. 

The last statistic to be discussed is the composite 
PDF for the first-level wavelet transform (Fig. 9). The 
wavelet transform is sensitive to inhomogeneity at a 
particular scale, and therefore, the PDF of the wavelet 
transform is related to the distribution of inhomoge- 
neity at the corresponding scale. Using a 1-h analysis 
window, the composite PDF (Fig. 9a) is symmetric 
about the mean and has a shape similar to a normal 
distribution. Since the PDF is peaked near the mean, 
this distribution has positive kurtosis, which increases 
with window size. We interpret this to mean that the 
contribution of inhomogeneity at a particular timescale 
is not evenly distributed over the window but is con- 
centrated to smaller regions that contribute a dispro- 
portionate amount of inhomogeneity to the window. 
This description corresponds to the concept of an in- 
termittent geophysical signal. As the window size in- 
creases to 8 h, the localized contributions of inhomo- 
geneity are better sampled and the PDF becomes more 
kurtose (Fig. 9b). An 8-h window violates the station- 
arity assumption and therefore raises the question 
whether the increase in kurtosis is an indication of 
cloud development or of better sampled intermittency. 

Results from the microwave radiometer LWP anal- 
ysis appear to be reasonable and provide a foundation 
for comparison to inhomogeneity model results. Both 
the ratio of standard deviation between scales and the 
slope of the power spectrum are time independent and 
appear to be related to each other, although the rela- 
tionship has yet to be developed. The correlation of 
standard deviation between adjacent levels is strong 
and provides motivation for the development of a 
wavelet-based cloud inhomogeneity model. This wave- 
let-based model must also include the Markov param- 
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(a) 1-h and (b) 8-h window sizes. 
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Fig. 9. LWP composite PDFs for first level detail with (a) 1-h and (b) 8-h window sizes. 


eter since it indicates the presence of correlated struc- 
tures at all timescales. On the smallest scale the 
distribution of the LWP inhomogeneity is similar to a 
normal distribution. Only as the window size includes 
a significant amount of cloud development does a 
strong positive skewness appear. This provides a basis 
to determine the suitability of an inhomogeneity model 
for a particular application. 

4. Inhomogeneity models 

The LWP inhomogeneity models used in this paper 
are multiscale statistical models. Operating within the 
framework of multiresolution analysis, values for the 
detail information at each scale are statistically gener- 
ated and then combined using an inverse wavelet trans- 
form. The detail information d m (t n/ 2 >») at scale m is 
dependent on the smoothed information s >,((/ j 2 m ) of the 
same scale and the detail information of lower resolu- 
tions d m + l (t n/ 2 ™ + >), d m+2 (t„, . . .. The values used 
to initialize the lowest resolution and the relationship 
between each scale provide the degrees of freedom for 
fitting the model. Cahalan’s BC model uses the values 
in s m (t n/2 >") to generate d m (t„, 2 »), while the VC model 
uses d m+ 1 ( t„ / 2 "’ + ' ) to generate d m ( t„ , 2 »> ) - Although each 
model uses a different approach, both still fit within the 
multiscale statistical framework. 

a. Bounded cascade model 

The BC model is a two-parameter model that uses 
an iterative cascading process to generate inhomoge- 
neity. Although the model is more thoroughly de- 
scribed in Cahalan et al. ( 1994a), a brief discussion is 
provided here. The BC model begins as a uniform 
cloud of thickness W, which is equal to the mean of 
the LWP. The first step of the cascade model moves a 
fraction of LWP,/ 0 , from one side of the cloud to the 
other. The direction that the LWP fraction is moved is 
arbitrary. Upon completion of the first level of the cas- 
cade the cloud consists of two regions. The first has a 
uniform thickness of W( 1 - / 0 ), while the second has 
a uniform thickness of W ( 1 + / 0 ) . The cascade process 
is repeated by moving a fraction /, from one side to the 
other for each of the two new regions of cloud. This 


results in four regions that have the following liquid 
water amounts: VL(1 +/ 0 )(1 +/i),W(l +/ 0 )( 1 — /), 
W(l -/ 0 )(1 + / ) , and W ( 1 -/>)( 1 -/,). Each time 
this process is repeated the number of regions is dou- 
bled, thus adding inhomogeneity to the cloud. Since the 
LWP is being redistributed, the mean LWP is con- 
served, and only the distribution of LWP is changed at 
each cascade. If the fraction moved at each cascade is 
the same, f 0 = /, = f 2 = • • • , after an infinite number 
of cascades the model will consist of an infinite number 
of singularities with varying intensities. Although this 
singularity model is truly a self-similar model, it gen- 
erates results that are nonphysical. This model is made 
more physically sound by adding a second parameter 
c that reduces the fraction moved at each level of the 
cascade. In this case 

fa = fc°; f =fc‘; f f,=fc" 

where c < 1. (8) 

No singularities are generated for this model since an 
insignificant amount of LWP is redistributed after a 
large number of cascades. By adding this damping fac- 
tor the model is no longer self-similar but multiaffine 
since there is a rescaling of the amount of LWP moved 
at each step of the cascade (Marshak et al. 1994) . This 
damping factor directly affects the amount of variabil- 
ity at each scale and therefore affects the slope of the 
power spectrum. The PDF of LWP generated by the 
BC model is lognormal in shape. 

The multiscale relationship of the BC model is illus- 
trated in Fig. 10 using the multiresolution framework. 
The initial uniform cloud LWP, W, corresponds to 
2<), which is the lowest resolution of the model. 
The detail d,-(?, i/2 .) is generated by multiplying 
Sj(t n , 2 0 with/o, which is the fraction moved by the first 
cascade. The smoothed and detail information are com- 
bined by performing an inverse transform with the Haar 
wavelet. This process is equivalent to generating a sin- 
gle level of the cascade model. The detail information 
for the second cascade is generated by multiplying val- 
ues of the smoothed information 57 (t n , 2 >- 1 ) with/, and 

repeating the inverse transform. Any level of inhomo- 
geneity can be generated by repeating this process 
through the desired number of iterations. The sign of 
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Fig. 10. Diagram of the BC model cast into the multiresolution framework. 


each coefficient in the detail information is randomly 
generated and determines whether the liquid water frac- 
tion is moved from right to left (positive) or from left 
to right (negative). This model is classified as a mul- 
tiscale statistical model on the basis of this random el- 
ement. Although the PDF, spectra, and moments of the 
model are deterministic, each realization of the model 
is randomly generated. 

Within the multiresolution framework, a number of 
variants of the BC model become apparent. The sim- 
plest variant comes with the choice of wavelet used in 
the inverse transform step of the cascade. Since the BC 
model redistributes LWP based on a Haar wavelet 
function, discontinuities between adjacent regions are 
generated as seen in Fig. 11a. By using a wavelet with 
a higher-order approximation, such as the D8 wavelet, 
the discontinuities at the boundary between adjacent 
regions are removed and the smoothness of the resul- 
tant LWP values is increased (Fig. 1 lb) . As a result of 
this change, the maximum value is reduced and the 
skewness of the generated LWP is also reduced. A dis- 
advantage of this change is the presence of negative 
LWP. This represents a nonphysical condition that 
must be modified in order to apply the results to radi- 
ative transfer problems. 

b. Variance coupled model 

Two additional changes are made to the BC model 
in order to define the VC model. Unlike the BC 
model, where the inhomogeneity at one level is de- 
pendent on the local mean, the VC model uses the 
inhomogeneity at the previous level to generate the 
detail information at the next level. As indicated in 


the analysis of the LWP data in section 3b, the stan- 
dard deviation at any particular level is well corre- 
lated with the standard deviation at the previous 
level. Although the mean is related to the standard 
deviation at a particular level, it is not as well cor- 
related and has a time-dependent element. When in- 
itializing the VC model with a uniform cloud layer, 
the mean LWP is used to generate inhomogeneity at 
the lowest level. Once this is done, the remaining 
levels of inhomogeneity are generated using only the 
inhomogeneity at the previous level. 

The second change in the BC model is related to the 
method of redistributing the LWP. The BC model re- 
distributes LWP using a deterministic fraction of the 
local LWP. The random element of the model is a bi- 
nomial value of either - 1 or + 1, which determines the 
direction that the LWP is moved. The VC model uses 
a Gaussian distribution of mean zero to redistribute the 
LWP. The standard deviation of this distribution is de- 
termined by using the ratio of standard deviation be- 
tween levels and the standard deviation of the previous 
level. 

The VC model was developed using multiresolution 
analysis and is best presented within this framework. 
The model is initialized by prescribing a series of val- 
ues to Sj(t n/2 i)- As illustrated in Fig. 12, the detail at 
the lowest resolution is generated from this initializa- 
tion data. The detail is assumed to have a Gaussian 
distribution of mean p and standard deviation a. Since 
the mean of the wavelet transform is near zero, only 
the standard deviation needs to be determined. The 
standard deviation of the detail d^t n , 2 <) is not constant 
and is calculated by multiplying a coupling factor 
with the corresponding value in Si(t n/2 ‘). The detail at 
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Fig. 11. LWP field generated from the BC model using the (a) D2 (Haar) and (b) D8 wavelets. 


the lowest resolution is then generated using the equa- 
tion 


di(t n ) = G(0, + a-diit „_,), (9) 

where G(/x, a) is a random number with Gaussian dis- 
tribution and a is an autocorrelation factor. The auto- 
correlation term is included because wavelet analysis 
of the LWP data indicates that a Markov process is 
present at each level of the wavelet transform. 

To generate the variability at the next highest level, 
the standard deviation is calculated using the variability 
of the current level. The standard deviation of the cur- 
rent level is multiplied by the variance coupling param- 
eter C a , which is the coupling between standard devi- 
ations of adjacent levels of the wavelet transform. The 
detail is then generated with the equation 

dj — \[t2n+k) G(0, C a ‘ (T j ( t n ) ) T OL dj \ ( ?2«+A:— 1 ) 

(* = 0 , 1 ), ( 10 ) 


where cr, (0 is defined as the standard deviation of a 
window of width 2k centered about t„ in the series 
di(t„); at {t n ) is calculated using the equations 

<r?(<„) = T- 1 — x (4U, + ;) - MO) 2 (11) 

2a - 1 
and 


1 V 

di ( L ) ^ dj(t„+j), 

A j=-x 


( 12 ) 


where ^,(J„) is the mean over this same window. It is 
necessary to calculate the standard deviation over a 
window of d, (r„) in order to provide a sufficient sam- 
pling and yet allow for local fluctuations in the standard 
deviation. Additional levels of detail are generated us- 
ing Eq. (10) since the coupling with respect to standard 
deviation appears to be independent of scale. As in 
multiresolution analysis, the detail information at each 
level is combined with the initialization data to generate 
the full resolution of the model. The choice of wavelet 
for the inverse transform is an additional variable. 
However, to be consistent with the analysis of the LWP 
data, the D8 wavelet is used. 


c. Analysis 

Analysis of the BC and VC models corresponds with 
that of the LWP data. The primary parameters used to 
generate the models’ results include the mean, standard 
deviation, and slope of the power spectrum. Values for 
these parameters are recorded in Table 2 along with 
other relevant parameters. An autocorrelation parame- 
ter is also included for the VC model since this is not 
inherently generated as in the cascade process of the 
BC model. 

Performing a windowed analysis on the BC model 
demonstrates the correlation between mean thickness 
and standard deviation. As seen in Figs. 13a,b, the 
mean and standard deviation appear identical apart 
from a scaling factor. This is an inherent property of 
the BC model where the variance parameter / deter- 
mines the ratio between these two statistics. Since in- 
homogeneity is generated at each scale of the model 
through the use of the local mean, correlation of stan- 
dard deviation between scales is also quite strong. The 
ratio of standard deviation between levels 1 and 2 is 
fairly uniform apart from a certain amount of noise. A 
similar relationship exists between standard deviations 
at larger scales. Although the inhomogeneity at any one 
level is a nonlinear combination of the mean and in- 
homogeneity at all previous scales, the ratio of standard 
deviation between levels is strongly influenced by the 
spectral parameter. 

Although the analysis of the VC model appears 
markedly different from the BC model, the relationship 
between scales is still present. The most prevalent dif- 
ference between the two models is that the uniform 
initialization field of LWP in the VC model is not 
greatly redistributed (Fig. 13c). This is due primarily 
to the model being initialized with 64 data points each 
having a value of 100 g m 2 . The standard deviation 
(Fig. 13d) does show variability; however, the mag- 
nitude is not as great as for the BC model. Since this 
model does not force a redistribution of LWP, the fluc- 
tuations in standard deviation give an indication of the 
natural variability of standard deviation within the VC 
model. 

The standard deviation at each scale of the VC model 
correlates well with the standard deviation at other 
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scales. This is a natural outcome of the VC model’s use 
of the variance coupling parameter. Although the pa- 
rameter is set at a value of 0.9, the calculated ratio of 
variability between scales is 0.99. This indicates that 
the ratio of standard deviation between scales is not the 
best estimator of the variance coupling parameter and 
a better estimator is needed. One advantage of the vari- 
ance coupling parameter is that the ratio of standard 
deviation between each level remains constant. This 
does not exist in the BC model since the ratio decreases 
as you go to higher levels in the analysis. This differ- 
ence exists because the D8 wavelet used in the VC 
model generates a much smoother field than the D2 
wavelet of the BC model. The discontinuities generated 
in the BC model introduce higher frequencies at each 
level of the cascade process. This artificially enhances 
the high-frequency component, and therefore, the ratio 


of standard deviation between scales approaches unity 
at small scales. With the VC model there is not as much 
frequency leakage and the ratio of standard deviation 
is maintained to the smallest scales. 

The power spectrum of both models is time inde- 
pendent and its slope is determined by the spectral pa- 
rameter of each model. For the spectral parameters used 
in Table 2 the slope of the power spectrum is —1.7 and 
- 1.6 for the BC and VC models, respectively. The dif- 
ference in the values used for the spectral parameter of 
the two models is related to the method of generating 
the inhomogeneity and the wavelet used in the inverse 
transform. The choice of wavelet has the most profound 
effect since it is not possible for a Haar wavelet to 
generate a power spectrum slope more negative than 
—2 since it is a step function (Barker and Davies 
1992). For particular ranges of the spectral parameter. 


Table 2. Parameters for the BC and VC models. The number of parameters specified gives an indication 
of the flexibility of the multiresolution framework for modeling. 


Quantity 

Bounded cascade 

Variance coupled 

Mean 

W = 100 g m -2 

sjx) = 100 g m 2 

Variance parameter 

/= 0.6 

C„ = 0.06 

Spectral parameter 

OO 

o 

11 

C„ = 0.9 

Autocorrelation 

NA 

a = —0.3 

Wavelet type 

Haar (D2) 

D8 

Coupling window width 

1 

\ = 8 

Number of points at 



initialization 

1 

64 

Number of cascades or 



levels 

i = 14 

i = 8 

Redistribution statistic 

Binomial -* either — 1 or + 1 

Gaussian-mean of zero 
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Fig. 13. Position dependence of the BC model’s local (a) mean and (b) standard deviation and the VC 
model’s local (c) mean and (d) standard deviation. 


the VC model as well as a modified BC model using 
higher-order Daubechies wavelets are able to generate 
power spectral slopes more negative than —2. 

As with the slope of the power spectrum, the auto- 
correlation of the wavelet transform at each level is 
time independent. Results from each model have an 
autocorrelation of -0.31. This corresponds to the value 
of the Markov parameter used in the VC model. The 
BC model generates autocorrelation through the cas- 
cading process and therefore has an advantage over the 
VC model since an additional parameter is not needed 
to determine the autocorrelation. 

The skewness of each model poses the greatest differ- 
ence between the BC and VC models. As mentioned pre- 
viously, the BC model inherently generates skewness 
through the cascading process. Using the parameters in 
Table 2, the skewness is 0.45, although it can be de- 
creased by decreasing the values of both the variance and 
spectral parameters. For the VC model the skewness re- 
mains near zero. Apart from conditions that generate neg- 
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Fig. 14. Composite PDFs of the 


ative LWP, the skewness of the VC model is unaffected 
by the choice of the variance and spectral parameters. 

The effect of skewness is easily seen in the composite 
PDF of results from both models. As illustrated in Fig. 
14a, the composite PDF for the BC model is clearly non- 
Gaussian and is actually lognormal as a result of the cas- 
cading process. This shape is present at all scales since the 
fundamental cascading process does not change for smaller 
scales. On the other hand, the composite PDF for the VC 
model is symmetric about the mean and is nearly Gaussian 
for all temporal scales (Fig. 14b). This demonstrates a 
clear difference between the two models and is partially 
due to the choice of wavelet and is mostly dependent on 
the redistribution process of the model. 

5. Discussion and summary 

a. Liquid water path data 

The analysis of the LWP data reveals some relation- 
ships that are useful for evaluating the effectiveness of 
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from the (a) BC and (b) VC models. 
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the two models. The most important is the relationship 
between inhomogeneity at different scales. From the 
windowed analysis it is evident that the inhomogeneity, 
quantified by the standard deviation, is correlated to the 
local mean of the LWP data. This inhomogeneity is 
distributed over many scales and, through the use of 
wavelet analysis, is seen to have a set relationship be- 
tween scales. This relationship is quantified by defining 
the ratio of standard deviation between scales. Since 
this ratio has a value that is less than 1 , the amount of 
inhomogeneity decreases at smaller scales. This ratio 
is present in the LWP data as well as in the BC and VC 
models, and indicates the existence of scaling proper- 
ties. Being time and scale independent, the ratio of stan- 
dard deviation between scales is closely related to the 
slope of the power spectrum. 

Another important relationship is the presence of 
skewness in the LWP data and its dependence on win- 
dow size. Using a half-hour window (9-km spatial 
scale), the skewness is very small and the composite 
PDF is nearly a Gaussian distribution. This Gaussian 
shape may not be indicative of the LWP data due to 
instrument noise and limitations in the retrieval algo- 
rithm. However, limiting measurements to those where 
the expected noise is less than 25% of the variability 
does not appreciably change the PDF. This implies that 
either the inherent skewness at the small scale is rela- 
tively small or the noise in the data is underestimated. 
With the 2-h window (36-km spatial scale) the skew- 
ness increases and noticeably affects the composite 
PDF. One can only speculate on the source of this 
skewness; however, a few observations can be made 
when comparing the LWP data with model results. Un- 
like the BC model, the presence of skewness is not 
present at all scales. This implies that certain dynamic 
processes that introduce skewness become important 
only beyond the half-hour timescale (9-km length 
scale). By comparing the time series for skewness and 
the mean, the skewness appears to be largest when 
cloud thickness is changing. This may indicate that the 
thickening and thinning process of cloud dynamics is 
nonlinear in marine stratocumulus. 

b. Inhomogeneity models 

Evaluation of the BC model shows that the model is 
useful in capturing some very important features of the 
data. This model is able to capture the correlation be- 
tween standard deviation and mean as well as the cor- 
relation of standard deviation between scales. Using 
only two parameters, the BC model not only adds in- 
homogeneity consistent with the standard deviation and 
power spectrum but also reproduces the autocorrelation 
and large-scale skewness that leads to a lognormal dis- 
tribution of LWP. The limitations of this model are the 
discontinuities generated by using a step function in the 
cascade process and the presence of skewness to the 
smallest scales of the model, which appears to conflict 
with results from the LWP analysis. 


Results from the VC model demonstrate both 
strengths and weaknesses of this model. As with the 
BC model, the VC model is able to capture the corre- 
lation between variability at different scales. However, 
due to the redistribution process of the VC model, the 
standard deviation is not as closely correlated to the 
mean. The standard deviation, power spectrum, and au- 
tocorrelation can all be adjusted using the parameters 
of the model. The VC model is best suited for gener- 
ating small-scale inhomogeneity within a known large- 
scale cloud field. Although this model produces distri- 
butions that are nonskewed locally, it is not able to 
generate the skewness that is present at the larger scales 
of the LWP data. A possibility for future work is to 
make a time-dependent version of the VC model that 
is able to generate skewness through the use of a non- 
linear cloud development scheme. 

The difference between the BC and VC models dem- 
onstrates the flexibility of the multiresolution frame- 
work for modeling cloud inhomogeneity. A number of 
variants to the BC model are implied by this research, 
such as the change in the cascading wavelet and the 
random number generator function. Although a high 
correlation' exists between inhomogeneity at adjacent 
scales, the method for generating this correlation can 
take a variety of forms. Finally, the multiresolution 
framework has an analogous two-dimensional form 
that provides a means of extending these inhomoge- 
neity models into higher dimensions. 

Both models involved in this study use fractal prin- 
ciples to generate cloud inhomogeneity. As seen in the 
analysis of the LWP data and the inhomogeneity mod- 
els, one way in which this principle manifests itself is 
in the correlation of inhomogeneity at different scales. 
This provides not only a method of generating a con- 
sistent PDF but also a constraint when adding inho- 
mogeneity into a known cloud LWP field. This feature 
may not be important in general circulation models 
where knowledge of the PDF appears to be sufficient 
for calculating radiative fluxes (Cahalan 1994b) . How- 
ever, consistent spatial distribution becomes important 
when calculating the radiance of a cloud field for re- 
mote sensing applications. In this later application the 
BC and VC models may prove useful for incorporating 
subresolution cloud inhomogeneity. 
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